Coarse-grained computations of demixing in dense gas-fluidized beds 
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We use an "equation-free", coarse-grained computational approach to accelerate molecular 
dynamics-based computations of demixing (segregation) of dissimilar particles subject to an up- 
ward gas flow (gas-fluidized beds). We explore the coarse-grained dynamics of these phenomena 
in gently fiuidized beds of solid mixtures of different densities, typically a slow process for which 
reasonable continuum models are currently unavailable. 
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Particulate flows, even for experimental systems of 
small size (~ 10 cm), consist of a very large number of 
discrete dissipative particles. Molecular dynamics (MD) 
simulations often serve as a quantitative modeling tool 
for such flows; however, such simulations for realistically 
large temporal and/or spatial scale problems are chal- 
lenging even with modern computers. Navier-Stokes-like, 
macroscopic continuum models have been developed by 
many authors based on the kinetic theory of granular ma- 
terials (see Refs. [H, [|| and references therein); however, 
quantitative continuum models for realistic particles (ac- 
counting for frictional interactions, heterogeneity among 
particles, and/or other inter-particle forces, such as van 
der Waals forces) in many regimes of practical interest 
(e.g. dense and/or cohesive flows where enduring con- 
tacts between particles occur) are currently unavailable. 

In this paper, we consider well-known phenomena for 
which the derivation of continuum models is still in flux; 
mixing and demixing (segregation) can occur when dis- 
similar particle mixtures of different sizes and/or densi- 
ties are subject to a strong enough upward fluid flow Q. 
A few different continuum models, more phcnomcnolog- 
ical or more rigorous, have been proposed 0, Hj], which 
often reproduce the phenomena in a qualitatively cor- 
rect manner' however, quantitative agreement is gener- 
ally elusive [5[. Furthermore, kinetic theory-based con- 
tinuum models for binary mixtures are much more com- 
plicated than those for uniform particles, and numerical 
simulation becomes more time-consuming (e.g. by an 
order of magnitude in Ref. ||). Accelerating the com- 
putation using (quantitative) microscopic models would 
therefore be invaluable for such problems. The objec- 
tive of this paper is to demonstrate a multi-scale com- 
putational approach enabling accelerated integration of 
MD-based microscopic simulations of dense particulate 
flows. 

Model. The particles are modeled as uniform-sized soft 
spheres (which can have different mass), whose inter- 
particle contact force F cont is determined following a 
model of Cundall and Strack @. The gas phase hydro- 



dynamics is accounted for in a volume-averaged way 0] • 
This approach has been used to study size difference- 
driven demixing [f|; here we consider density-driven 
demixing. 

We deliberately choose demixing occurring in narrow 
beds (cross sectional area of 15<i p x 15d p with periodic 
boundary conditions for both lateral directions, where d p 
is the particle diameter) as a test problem, so that the ex- 
act results can readily be computed and used to critically 
test the coarse-grained computations. These are quasi- 
1D flows, where the coarse-grained gas flow is effectively 
ID, while particle simulation is fully 3D. Demixing be- 
comes more pronounced Q in such narrow beds. The 
equation of motion for each particle can be simplified to 
be O: 
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where m p and v p are individual particle mass and veloc- 
ity, respectively; g is the gravitational acceleration; V p 
is the volume of each particle; <j> is the locally-averaged 
solid phase volume fraction; J3 is the interphase momen- 
tum transfer coefficient [ll|> [HI ; Us is the coarse-grained 
solid phase velocity; and U s is the superficial gas flow 
velocity, in the direction opposite to that of the gravity. 
In our study, the Reynolds number based on the particle 
size is generally very small (<~ 0.1), and (3 is approxi- 
mated to be 
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where fi g is the gas phase viscosity. We nondimensional- 
ize quantities by using p 



dp, \J gd p and \fdp/g as char- 
acteristic density, length, velocity, and time scales, where 
p s is the solid phase mass density of the lighter particles. 
More details of the model can be found in Refs. (Tol 

Direct simulations. Demixing is typically a slow pro- 
cess, whose occurrence and duration depend on the den- 
sity difference and the gas flow rate. Direct simulation 
with a sufficiently large U s (= |U S |), starting from a ho- 
mogeneously mixed, packed (static) state (Fig. Q]), illus- 
trates that particles of different densities gradually demix 



FIG. 1: (Color online) Snapshots of a gas-fluidized bed of a 
binary mixture of identical size but different density particles, 
undergoing spontaneous demixing, shown at times equally 
separated by At = 100. Light-colored (yellow) particles are 
twice as dense as the dark-colored (red) ones (coefficient of 
restitution (friction) = 0.9 (0.1); k n = 2.0 x 10 5 ; U, = 0.41; 
d p = 100 /mi; p s = 0.90 g/cm 3 ). 



spontaneously. When U s is well above the minimum 
fiuidization rate of both species (as in Fig. [I}, the bed 
exhibits ID traveling waves (TWs) [13j, and demixing 
occurs superposed on the persistent oscillatory motion 
driven by lD-TWs. A typical computation of an entire 
demixing process in the above "tiny" system (2 x 10 7 in- 
tegration steps of 12 500 particles shown in Fig. []} takes 
nearly two days, or more than a week for smaller gas 
flow rates, on a single-processor PC with 1.7 GHz CPU. 
Obtaining an ensemble of long simulations for statistical 
averaging purposes can be extremely time-consuming. 

Coarse-grained description and "observables". In the 
literature, the degree of mixing/demixing is often char- 
acterized by various lumped indices (or "order parame- 
ters" ) pj . [14 HEl HU such as the so-called Lacey mixing 
index [lj](e.g. see its use in Ref. 8j). Here we seek 
coarse-grained variables (or "observables" ) that could be 
used in continuum descriptions. 

As in the two-fluid modeling approach [l| , it is natural 
to think of hydrodynamic variables as candidate coarse 
observables. From direct simulations, we observe that 
in the course of demixing in quasi-lD beds, the process 
strongly depends on the local density, which makes the 
ID volume fraction profiles themselves sufficient coarse 
observables; when we suddenly randomize only the in- 
dividual particle velocities (hence the granular tempera- 
ture as well) and continue the simulation, the demixing 
progresses essentially undisturbed. We further recognize 
that cumulative particle distribution functions (CDFs, 
and more precisely, their inverses which are bounded by 
and 1) are more convenient coarse observables: CDFs 
are smoother than volume fraction profiles, suffer from 
less noise, and facilitate the lifting (see below) procedure. 

We consider both discretized inverse CDFs (ICDFs; 
Fig. [2]) and their compact parametric representations 
(finite element or expansion coefficients in convenient 
polynomial sets; Fig. [3} as our actual coarse observ- 
ables. For a narrow range of U s 's slightly above the 
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FIG. 2: (Color online) Direct, full integration (thin lines) 
and coarse projective integration (patches of thick lines), us- 
ing evenly discretized CDFs as coarse observables. (a) Slow 
demixing at Us = 0.19, where the bed hardly expands and 
no lD-TWs form. Projective steps of At = 80 were taken 
through forward Euler, after direct integration for 8t = 20 
(the latter half of data were used to estimate the local slope), 
(b) In the presence of lD-TWs (U s = 0.41); 8t = 30, and 
At = 2T, where T is the average period estimated during 
short bursts of direct integrations; the last two periods of the 
locally oscillating data were used to estimate the coarse slope. 



minimum fiuidization rate, the time evolution of ICDFs 
does not exhibit any waviness (thin lines in Fig. [2] (a)), 
and demixing occurs very slowly (complete demixing 
is hardly achieved). In the presence of lD-TWs at 
larger J7 s 's, ICDFs locally oscillate at a fast time scale 
(Fig. [21 (b)); these oscillations can be smoothed through 
ensemble-averaging. Once the coarse observables are cho- 
sen, governing equations for their time evolution need to 
be derived for further computation. We will follow an 
equation-free approach, circumventing the derivation of 
such equations, assuming they do exist conceptually, but 
are not available explicitly. 

An equation-free approach. When the time series of 
coarse observables (obtained by direct integration of the 
microscopic simulator) are smooth and slowly varying, 
one can estimate their local time derivatives and then 
project the values at a future time (e.g. using forward 
Euler or more sophisticated schemes) . We recognize that 
if one can initialize the microscopic simulator consistent 
with the future (projected) values of the coarse observ- 
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ables, one can actually accelerate the overall computa- 
tion. This simple idea underpins coarse projective inte- 
gration [13, Cl- 
in equation-free computations, traditional continuum 
numerical techniques are directly applied to the out- 
come of appropriately initialized short bursts of micro- 
scopic simulation, and the macroscopic equations are 
"integrated" or "solved" without ever being written 
down [lij [111 The essential steps are: (i) Identify 
coarse observables (which are discretized ICDFs or their 
parameterization coefficients in our study). For conve- 
nience, we denote the microscopic description (here the 
individual particle positions) by x, and the macroscopic 
description (here the ICDFs) by X. (ii) Choose an ap- 
propriate lifting operator /i£ , which maps X (ICDFs) to 
one (or more, for the purposes of variance reduction and 
ensemble-averaging,) consistent description(s) x (here, 
particle positions). Figuring out an efficient lifting op- 
erator is essential, (iii) Starting from lifted initial condi- 
tion^) x.(to) — /xi,(X(fo)), run the detailed simulator for 
some time horizon (X^ > 0) to obtain x(4q + Tft). (iv) 
Use an appropriate restriction operator A4r which maps 
the microscopic state(s) to the macroscopic description 
X(t + T h ) = Mn(x(to + T h )), resulting in time series of 
the coarse observables (ICDFs), or a coarse time-stepper 
$ Th for them: X(t + T h ) = $ Th (X(t )). (v) Apply de- 
sired numerical techniques (forward Euler, in our study) 
to the coarsely observed results in step (iv), and repeat. 

Lifting. Given an ICDF as the coarse observable, we 
need to construct particle configurations consistent with 
it. Arranging particles (i.e. sphere packing) in three- 
dimensional space with an arbitrarily prescribed ICDF 
(or local volume fraction, especially dense) profile is non- 
trivial and generally time-consuming [22j |. as excessive 
particle-particle overlap has to be avoided; it would be 
less difficult for dilute particulate flows. In our study, the 
total height of a bed remains virtually the same, even in 
the presence of 1 D-TWs (see the second and later frames 
in Fig. [T]). Therefore, we do not reassign particle loca- 
tions from scratch in each lifting step. Instead, we uti- 
lize particle locations which are already obtained from an 
earlier step, and switch only particle indices (or "colors" , 
whether "red" or "yellow" ) until the prescribed ICDF be- 
comes satisfied. When a system of polydisperse particles 
is considered, this scheme should be modified; this would 
be a subject of future research. 

Our lifting operator for single realization computa- 
tions requires additional consideration when ID-TWs are 
present: particle locations obtained from earlier simula- 
tion at the same "phase angle" during the wave propa- 
gation have to be used. In ensemble-averaged computa- 
tions, the oscillatory motion becomes smoothed, requir- 
ing no special attention. 

Coarse projective integration. We choose discretized 
ICDFs of the lighter particles as the coarse observables, 
and accelerate the demixing computation using coarse 
projective forward Euler scheme [13, EH ■ For smaller U s , 
where the bed hardly expands and ICDFs do not oscil- 
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FIG. 3: (Color online) Evolving inverse CDFs (ICDFs) for 
the same case as in Fig. [T] averaged over 10 realizations: (A) 
ICDF of particles physically located in the upper half of the 
bed, and (B) ICDF of the lighter particles, which are shown at 
t = 150. Solid lines with negative slopes are snapshots of the 
difference between (A) and (B) at different times; compare 
with the functional fit shown as dashed lines; see Eq. ©• 

late, the projection step size is determined by only the 
temporal smoothness (accuracy of the local linearization) 
of ICDF evolution. The demixing occurs very slowly in 
this case (Fig. [5] (a)), and these computations can achieve 
high computational speedup. Excessively large projec- 
tion steps can cause inaccuracies, similar to large time 
steps in normal integration. In the presence of ID-TWs, 
the projection step size is chosen to be an integral mul- 
tiple of the local oscillation period (Fig. [2] (b)). 

Projectively integrated values (thick lines in Fig. [2]) 
follow the trajectories of direct, full integrations (thin 
lines in Fig. \2§ well, confirming that these coarse observ- 
ables are good continuum variables. Ensemble- averaging 
of ICDFs over different realizations (of different phase 
angles during wave propagation) smoothens local oscil- 
lations arising from ID-TWs. Projective integration of 
ensemble-averaged ICDFs, both in the presence and ab- 
sence of ID-TWs, can be applied in the same way. 

A more compact description. The difference between 
the ICDF of particles located in the upper half 'of the bed 
(irrespective of their densities) and that of the lighter 
ones can serve as a useful coarse observable; ICDFs of 
the lighter particles are bounded by those of the (total) 
top half particles, and their difference is positive definite 
(Fig. [3]) . Furthermore, once nearly- full local demixing oc- 
curs, the ICDF difference there becomes virtually zero. 
The difference of the two ICDFs can be fit by the follow- 
ing simple functional form (dashed lines in Fig. [3]): 

y = max[Ax + B + Cexp(Dx),0], (3) 

where A(t), B(t), C(t), and D{t) are our new coarse ob- 
servables, determined on the fly by functional fitting, and 
x and y represent the abscissa and the ordinate in Fig. [3] 
respectively. Other basis functions, such as high-order 
polynomials also can fit the functional form reasonably 
well, but they require many more terms and the time evo- 
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FIG. 4: (Color online) Comparison between direct, full in- 
tegrations (solid lines) and coarse projective integrations 
(groups of dots), using the four coarse observables in Eq. ([3]). 

lution of their expansion coefficients is not generally slow. 
Lifting for these new coarse observables involves a minor 
intermediate step: mapping between these four variables 
and an ICDF discretization through the functional form 
in Eq. ©. 

We use these four coarse observables to perform 
ensemble-averaged coarse projective integration over a 
number of realizations (Fig. 2]) . These observables vary 
slowly and smoothly in time (occasional oscillations dis- 
appear for larger ensembles); in a sense, this is a pseu- 
dospectral solution of the unknown governing equations 
for ICDF evolution. 

Conclusions. We have used an "equation-free" coarse- 



grained approach to accelerate (by a factor of two to ten; 
the lifting step in our study involves minimal compu- 
tational effort) computations of dense particulate flows 
and, in particular, of demixing occurring in gas-fluidized 
beds of dissimilar particles. This approach holds promise 
for the prediction of coarse-grained behavior at practi- 
cally relevant spatial and temporal scales. 

We deliberately considered a quasi- ID illustrative 
problem in our study, in order to demonstrate the vi- 
ability of the approach. As a consequence of the problem 
considered in this study, the coarse observables were one- 
dimensional discretized ICDFs. For systems involving 
higher dimensional flows, candidates for coarse observ- 
ables may include marginal and conditional ICDFs [23|. 
More work for such systems has to be done to identify 
proper coarse observables and an efficient lifting opera- 
tor, vital components of this approach. Ensemble aver- 
aging reduces fluctuations among the realizations, giving 
better quantitative representations; the computation of 
each realization readily parallelizes across computational 
nodes. 

More sophisticated equation-free algorithms (e.g. 
coarse fixed point algorithms (2U^ ) can be used to find sta- 
ble as well as unstable steady states; quantify their stabil- 
ity; and perform numerical bifurcation analysis. Exploit- 
ing such tools to investigate the coarse-grained dynamics 
of mixing and demixing (and other particulate flow prob- 
lems) is the subject of current research. 

This research was partially supported by DOE, 
DARPA, and ACS-PRF. 
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